# Module membership and gene significance

library(WGCNA)

exprFilePath = "/Users/zyl/Workspace/hzau/MOIN-Extension-Github/WGCNA-Extension/FemaleLiver-Data/datExpr.csv"
traitFilePath = "/Users/zyl/Workspace/hzau/MOIN-Extension-Github/WGCNA-Extension/FemaleLiver-Data/datTraits.csv"
traitName <- "weight_g"
moduleName <- "brown"


datExpr = read.csv(exprFilePath, stringsAsFactors = FALSE, row.names = 1);
datTraits = read.csv(traitFilePath, stringsAsFactors = FALSE, row.names = 1);

nSamples = nrow(datExpr);
load("./network.RData")


# Define variable weight containing the weight column of datTrait
weight = as.data.frame(datTraits[[traitName]]);
names(weight) = "weight"

moduleColors = labels2colors(net$colors)
MEs0 = moduleEigengenes(datExpr, moduleColors)$eigengenes
MEs = orderMEs(MEs0)

modNames = substring(names(MEs), 3)
geneModuleMembership = as.data.frame(cor(datExpr, MEs, use = "p"));
MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples))

names(geneModuleMembership) = paste("MM", modNames, sep="");
names(MMPvalue) = paste("p.MM", modNames, sep="");
geneTraitSignificance = as.data.frame(cor(datExpr, weight, use = "p"));
GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples));
names(geneTraitSignificance) = paste("GS.", names(weight), sep="");
names(GSPvalue) = paste("p.GS.", names(weight), sep="");

column = match(moduleName, modNames);
moduleGenes = moduleColors==moduleName;

png(
  filename = paste("/Users/zyl/Workspace/hzau/MOIN-Extension-Github/WGCNA-Extension/FemaleLiver-Data/", traitName, "_", moduleName, ".png"),
  width = 13,
  height = 9,
  units = "in",
  res = 300
)

par(mfrow = c(1,1));

verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
                   abs(geneTraitSignificance[moduleGenes, 1]),
                   xlab = paste("Module Membership(", moduleName, ")"),
                   ylab = paste("Gene significance(", traitName, ")"),
                   main = paste("Module membership vs. gene significance\n"),
                   cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = moduleName)


dev.off()
